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I. INTRODUCTION 


As the world's fossil fuel resources are depleted, more 
emphasis is being placed on the breeder reactor as a poten- 
tial means of solving the coming energy crisis. While the 
development of new energy sources iS being pushed, equal ef- 
fort is being given to the maintenance of an environmental- 
ly clean world. To this end, the safety of breeder reactors 
1s receiving a considerable amount of attention before assum- 
ing that the breeder reactor is the answer to the energy 
Droblem. 

The Liquid Metal Fast Breeder Reactor (LMFBR) appears to 
be one of the most promising breeder reactors. Most engineers 
will concede there is little probability of a nuclear explo- 
Sion occurring in the operation of a nuclear reactor. Of 
major concern to engineers is the loss of coolant accident 
(LOCA) and the transient overpower accident (TOP). The pres- 
ent analysis is concerned with a TOP accident in a LMFBR. 

The analysis is formulated to model the dynamic response of 
the reactor fuel subassembly during the initial period of the 
postulated overpower transient. The primary consideration is 
given to the early response of this fuel subassembly to 
Various conditions of disturbances. The phenomenon which oc- 
curs after core disassembly (i.e., clad melting) is not the 
concern of this analysis. Only the time prior to clad melting 


is being considered. 
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No consideration is given here as to how the overpower 
transient occurs or to why the safety features of the reac- 
maeaid not operate property. I[t is postulated that the 
accident has occurred. In this analysis, the TOP accident 
is created by either a step increase in reactivity, a ramp 
increase in reactivity, or a combination of both. 

An inherent safety feature of most reactors, nuclear 
Doppler broadening feedback, is included in the dynamic model 
of the fuel subassembly. The Doppler feedback acts to reduce 
the effect of the excursion. Consideration of this feedback 
creates a non-linear system model which is described by a 
non-linear, initial-boundary-value problem. 

The conventional method of solution uses the standard 
point kinetics formulation. Recent studies have pointed out 
a non-negligible error in this model [1], particularly with 
asymmetric disturbances [2], or space-dependent feedback [3]. 
In Ref. [4], a somewhat novel approach of using the finite 
element method (FEM) for the space-time dependent solution 
of the reactor dynamics problem was demonstrated. The FEM 
is effective in handling these asymmetric disturbances and 
Space-dependent feedbacks. Therefore, the finite element 
method was used so that the spatial effects on the postu- 
lated problem may be studied further. 

The purpose of this work was to demonstrate further the 
applicability of the FEM to the non-linear reactor dynamics 
Problem as well as toinvestigate the dynamic response of the 


reactor fuel subassembly. The analysis required a novel 
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approach to handie the gap conductances present at the in- 
terfaces of the equivalent cell model of the fuel subassem- 


bly; this will be clarified in the analysis. 
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re OeSCRIPTION OF PROBLEM 


Eyeeeem?olCAL SYSTEM 

The typical Liquid Metal Fast Breeder Reactor (LMFBR) 
core consists of many hexagonal modules, each containing 
several hundred fuel pins. For this analysis, an equivalent 
cylindrical cell is used to model the fuel subassembly; see 
Figure 1. The use of equivalent cells as models for larger 
Systems has been common practice in nuclear analysis (i.e., 
the well known Wigner-Sietz method). In using an equivalent 
cell, the actual shape of the reactor core is not important, 
and theanalysis is applicable to any reactor which has the 
Same equivalent cell. 

The equivalent cell considered in this analysis, Figure 1, 
is fueled with enriched uranium dioxide, has a stainless 
Steel cladding, and has liquid sodium for a coolant. The 


dimensions used are 


fu 
ll 


Reese cil. 
b = 0.292 cm, 
O2365 “cm, 


C 


and H 


Soe <cMm, 

The gap between the fuel and cladding is very small and, in 
Fact, may be nonexistent as in bonded fuels. The dimension 
of this gap has been assumed negligible. The height, H, of 
the fuel rod is shorter than many proposed systems (Fast Flux 


Testing Facility and Clinch River Breeder Reactor). However, 
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to facilitate the numerical solution a smaller rod was used. 
The dynamic behavior prior to fuel pin failure for this sys- 


tem should be similar to the behavior of larger systems. 


The treatment of this problem in three dimensions would 
be prohibitive in computer usage. Therefore, azimuthal sym- 
metry is assumed, and the problem becomes a two-dimensional] 


ee lindrical (r,z) problem. 


pees 1 So TEM. MODEL 

The analysis considers the monoenergetic neutron diffu- 
Si0n approximation to model the transient neutron transport 
problem. A simple conduction-convection heat transfer model 
is used for the energy transport problem. The Pempenaenee 
in the model are directly coupled to the neutron flux through 
the nuclear heat generation within the fuel. The neutron 
population is in turn coupled to the temperature through any 
of a number of reactivity feedback mechanisms. The nuclear 
Doppler effect is perhaps the most important of these mecha- 
nisms since it provides a negative temperature coefficient 
which increases the inherent stability of the reactor. Prior 
to core disassembly and fuel melting, the nuclear Doppler 
effect is the most dominant feedback and is, therefore, the 
only feedback mechanism considered in this analysis. 

For irradiated, mixed-oxide fuels, a phenomenon of fuel 
restructuring has been commonly observed. This restructuring, 
essentially a change of phase of the fuel material, presents 
a unique heat transfer problem particularly during transient 


conditions. The problem has not been fully characterized and 
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is beyond the scope of this analysis. The fuel ree there- 
fore, assumed to be a homogeneous mixture of enriched uranium 
dioxide. 

At the fuel-cladding interface, there exists a gap which 
produces a thermal resistance. This thermal resistance is 
One of the most significant deterrents to the energy transfer 
to the coolant. The interface may be in physical contact or 
an actual gap may exist. The prediction of the thermal resis- 
tance is extremely complicated and must take into considera- 
tion many parameters: initial dimensions, type of bond, fill 
gas composition, fuel restructuring, fuel swelling, prior fuel 
life cycle, to mention a few. Reference [5] documents a com- 
puter program which attempts to predict the gap conductance, 
Naps This treats the thermal resistance at the interface in 
the same manner as a convection heat transfer coefficient when 
considering convection heat transfer. Since it is not the 
objective of this analysis to predict the gap coefficient, a 
representative set of values for gap coefficient, as given in 
Ref [6], is used in this analysis. These values are assumed 
to remain static during the transient. The gap conductance 
mite actually varies with time and will have an effect upon 
the transient, as noted in Ref. [7]. The prediction of this 
Variance was not considered important for this analysis; 
therefore, the static assumption was made. 

An average convection heat transfer coefficient was used 
to determine the heat transfer from the cladding to the cool- 
ant. The value used was determined from an empirical formula 


Given in Ref. [5] and repeated in Appendix C. 
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In this work, consideration is given to step and ramp in- 
creases in reactivity, although any reactivity transient may 
easily be considered. The step and ramp increases in reactiv- 
ity probably represent the most realistic physical reactivity 
inputs in a reactor. Once the reactivity has been inserted, 
the transient overpower excursion begins. Unless the Doppler 
feedback can override the inserted reactivity, the excursion 


will continue until there is physical core disassembly. 


om NUMERICAL SOLUTION 

The system of equations which models the proposed problem 
is a non-linear, initial-boundary-value problem. The conven- 
tional method of solution of the reactor dynamics is the point 
kinetics formulation. It was pointed out in Refs. [1], [2], 
[3], and [4], that there is a non-negligible error in this 
model, particularly under conditions of asymmetric disturbances 
or space-dependent feedback. Reference [4] demonstrates the 
Somewhat novel approach of using the finite element method 
(FEM) to solve the space-time dependent reactor dynamics prob- 
lem. As shown in Ref. [4], the FEM is quite effective in 
nandling localized perturbations and space-dependent feedback. 
In this work, only uniform disturbances were considered; how- 
ever, the feedback model was space-dependent. Therefore, the 
Finite element method is used to solve the non-linear, coupled, 
Space-time dependent neutronic and heat transport field equa- 
tions. The solution technique results in a large computer 
storage requirement; therefore, an optimum compact storage 
scheme, Ref. [8], is utilized for storage of the discretized 


Maurices . 
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Once the domain has been discretized by the FEM, the solu- 
tion of the resulting ordinary differential equations was to 


be accomplished by the implicit Gear's method, Ref. [9]. 
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Pe ae OUE DS De VeLORME NT 


A. NEUTRONIC ANALYSIS 

In this section the governing field equations for the 
neutron population (flux) for each of the three regions 
(fuel, clad, and coolant) of the domain will be formulated. 
The monoenergetic, diffusion theory will be used. 

Consider an arbitrary volume of material within a reac- 
tor. Applying the condition of conservation to the mono- 
energetic neutrons leads to the neutron equation of 
Seroinuity [10] 


antrot) = s(rjt) - r(r)o(r.t) - div d(r,t) (1) 


where 


= Sauda le pom nt 


t - time 
n(r,t) - neutron density 
S{r,t) - neutron production 
Zr) - neutron absorption cross section 
o(r,t) - neutron flux 
J(r,t) - neutron current 


The left-hand side of equation (1) represents the time rate 
of change of the neutron density which is related to flux, 


Ds. OY 





where 
v - neutron velocity 


On the right-hand side of equation (1), the first term is 
neutron production, the second term is a neutron loss through 
absorption, and the third term is a neutron loss through leak- 
age from the control volume. 

Using equation (2) and applying Fick's Law to the equation 


of continuity results in the classical neutron diffusion equa- 


tion 
ve(D(r)vo(r.t)) -Z,(r)o(rst)+s(r,t) = 2 3848) (3) 
where D(r) - neutron diffusion coefficient 


The vector notation used here is intended to include only two 
dimensions (r,z) since azimuthal symmetry has been assumed. 

Equation (3) is applicable to each of the three regions 
eeeresequivalent cylindrical cell. The subscripts F (fuel), 
c (cladding), and co (coolant), will be used to denote these 
regions. 

1. Fuel Region 

In applying equation (3) to the fuel region the 

material properties of the fuel must be used. Within the 
fuel the neutron souce term is due to the nuclear fission 
Process. During fission, neutrons are released as both 


Prompt neutrons and delayed neutrons so that 


S(r,t) = S (r,t) + Sp(rot) (4) 
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where S (r,t) - prompt neutron source 


Sy(r,t) -~ delayed neutron source 


The neutron sources are commonly represented as [10] 


and 
n 
Sp =e Ui; (6) 
T=] 
where 
kK. - infinite multiplication factor 
ra- (traction Of f1ISSIOM Neutrons whichedppear 


as delayed neutrons 


. - number of delayed neutron groups 
f 
} 
C. - concentration of delayed neutron precursors 
in the ith group 
AX. - decay constant of the delayed neutrons 


The space and time variables, r and t, will be dropped except 
where needed for clarification. 

The concentration of delayed neutron precursors, C., 
1S given by the following first order partial differential 


equation [10] 


aU. 
: 
opm So. Be ey oe (7) 
where 
B. - fraction of delayed neutrons which appear 


: as delayed neutrons in the ith group 


24 








The solution of equation (7) is 


ie 
= -rA.(t-t') i ! ° -r.t 
G. = Bf a7" 4 k(Pst')Dapopdt'+ Cle "1 (8) 
O 
C° - concentration of delayed neutron precursors 
Of the ith group at time zero 


Reference [10] develops an expression for the initial concen- 


tration of delayed neutrons 


Up ee ae (9) 


where 


oe = Ne laie finite, mui eipinica ten facto, 


»° - initial neutron flux 


Soueiming equations (3), (4), (5), (6), (8), and (9), yields 


the governing equation for the fuel region 


V-(D-Vo-) + Door lk (1-8) -1] 


n i 
+ DUA.[8, ey tees Mkeo(rst' )oet, edt! 
es 
dg | 
° “A.t5 _ I F 
a ae Pr e 1 ] = ¥ at t 1i0:) 


2. Cladding Region 


The cladding separates the fuel and coolant and con- 
tains the fuel and the fission by-products. Equation (3) 


Seovevihiceene meutron TluUx in the clad. Since the clad contains 


DNs 





mgmetissile material, the neutron source term is zero. The 


field equation is, then, 


_ | 
ee ete cee note (11) 


meee Coolant Region 


The annular region around the cladding of the equiva- 
lent cell contains the coolant. As in the clad, there is no 
fissile material in the coolant and, therefore, no neutron 
Source term. The equation governing the neutron flux in the 


coolant is 





ol Go 
a ear 72 aad E (12) 


Equations (10), (11), and (12) are the one-velocity, 
diffusion approximation used to model the neutron transport 
problem. 


feat inite Multiplication Factor 


The infinite multiplication factor, k may be ex- 


Oeeased as the infinite multiplication factor at time zero 
Meearteat transient), Kee, PlUS Thempostulated) redeting ty in 
Sertion (such as a step or a ramp), po, minus the change in 
the Doppler reactivity feedback, Adp- Other feedback mechan- 
isms are normally not as significant as the Doppler broaden- 


ing feedback prior to fuel melting and have been neglected 


in this analysis. Therefore, 


co 


pei 0 = 80, (13) 


Zc 


For a fast reactor, k° may be approximated as 


2 
ko =v git (14) 
a F 
where 
Vv - average number of neutrons released 


per fission 
Der - fission cross section of the fuel 


a - absorption cross section of the fuel 


The nuclear Doppler effect is a very important safety 
feature in a nuclear reactor. Nuclei in an atom are in con- 
tinual motion due to their own thermal energy. As a result 
of this motion, even when monoenergetic neutrons interact 
with the atom, there appears to be a spread in the energy of 
the neutron - the Doppler effect. It can be shown that the 
cross section of a resonance becomes less in magnitude and 
wider as the motion of the nuclei increases [10]. As the 
temperature increases, the motion increases, and the shape 
of a resonance cross section broadens. This broadening in- 
creases the average cross section, thus, providing a negative 
temperature coefficient. This effect is shown in Figure 2. 
It is this nuclear Doppler broadening effect which provides 
one of the few inherent, reliable, negative reactivity feed- 
backs which slows an overpower transient and possibly stops 


a mild overpower transient. 
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Figure 2. Doppler Broadening of a Resonance Peak 











The Doppler reactivity change with respect to fuel 


temperature changes should be written as [6] 


QO. 


p ; 
Seale 


a 


T 


= 


a ppl ee eye. (15) 


QO. 


where a, b, and c are parameters determined from experi- 
mental work and m is an integer. However, as noted in 
Ref. [6], a substantial amount of work has shown that 

ij “Po is very nearly constant over the temperature range 


under consideration. Therefore, the coefficients a and 


Cc have been set equal to zero, and b is defined as 
DS Np aa (16) 


The constant, Ky 1S commonly called the Doppler constant. 


Solving equation (16) for Py yields 


Pp = Oe oe ar le) 


where K is an integration which may be obtained from 


initial conditions. 


where 


° 


a Doppler effect at time zero 


Te - fuel temperature at time zero 


PLS, 





Substituting for K in equation (17) will give 


ob) - 92 = App = b In(T_/T2) (19) 


The infinite multiplication factor now becomes 


oot In(T-/Te) (20) 
The effect of delayed neutrons is small compared to 
the prompt neutron effect; therefore, it may be assumed that 
the Doppler effect on delayed neutrons is insignificant to 
the overall problem. The k, of equation (8) is, then, 
assumed to be k® . With this assumption and equation (20), 


the neutron diffusion equation in the fuel, equation (10), 


may be rewritten as 


V* (DeVoe) + Da poelk2(1-8)-1+0(1-8)-b(1-8) In(Tp/T2) J 


ra (ee evo m en tee a E 
om ece,2,, f e, Va [kZtelo-dt EC cle }= a 
i=] O 


(21) 


To facilitate the present analysis, the number of 
delayed neutron groups is taken as one averaged group. So 
that, A. and 6. become h and 8, respectively. This approxi- 
Mation should have little effect on the problem under con- 
Sideration. 

ee Loundary Condi tions 

The neutron diffusion problem involves solution of 

ene partial dafferential equations (11), (12), and (21). with 


the following boundary, interface, and initial conditions: 
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1) ae (605 Zsa) = 0 
2) op (a,z,t) = > (a,z,t) 
9p 9O 
3) De or, Kaz; t) = De a ( aujiziet ) 
4) 6, (bazst) = 05(b,z,t) 
ae 29 69 
5 ) ve ea | D scene) = Use are (Gieez 5 ¢ ) 
a 
6 ) sa (c,Z,t) = 0 
7) oe (rat Bet)= O.(rat Set) = oO o(rat got) = 0 
8) be (150) = o2 (rr) 
9) o. (7,0) = o2 (r) 
10) boo(he0) = o20(r) 


Boundary condition 1) results from the assumed azimuthal sym- 
metry. Interface conditions 2), 3), 4), and 5) are continu- 
meyeeconditions of the flux. Boundary condition 6) results 
from the use of an equivalent cell and basically indicates 
there is an equal number of neutrons transferred in and out 
of the cell at the outer boundary. This should be valid un- 
less the cell is located near the outer edge of the reactor. 
Boundary condition 7) is an assumption that the flux is zero 
at the axial boundaries of the cell. Initial conditions 8), 
Peed Gomare the assumed Initial distributions of the 


meltron Tf lux. 
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B. HEAT TRANSFER ANALYSIS 

In this section, the principle of conservation of energy 
will be used to formulate the governing field equations for 
the heat transport in each of the three regions. A simple 
heat conduction model] with convection heat transfer to the 
coolant is used to model the heat transport problem. A gap 
conductance mode] is used to describe the heat transport 
across the gap at the fuel-clad interface. 

1}. Fuel Region 

Conservation of energy within the fuel region yields 


the unsteady heat conduction equation with a generation term 


: : F 
Ve(ke(r)VTe(rt))#a(rot)=op(r)Coe(r) gee (net) (22) 
where 
ke (4) - thermal conductivity of the fuel 
Te (r,t) - fuel temperature 
Q(r,t) - nuclear energy generation per unit 
volume 3 


pe (1) - fuel density 


Syeh Ie) - fuel specific heat 


As in the neutronic analysis, the vector notation is intended 
to include only two dimensions, (r,z). The r and t will be 
dropped except where needed for clarification. 


The nuclear generation term may be expressed as 


q =e Ler o¢ (23) 


a 





where e - nuclear energy released per fission 


As can be seen, it is through the nuclear generation term 
that the temperature is directly coupled to the neutron 
flux. This coupling and the temperature dependent Doppler 
reactivity feedback combine to make the coupled problem non- 
linear. 

Substituting equation (23) into equation (22) yields 
the governing thermal equation for the fuel 


dT. 


Pr Sor ot ee) 


— 
oe 


Ve(keVT_) + e Dee oe 


gee Cladding Region 


Conservation of energy within the clad will yield 
the heat conduction equation. With nuclear generation, a 
relatively small amount (~5%) of the energy will be released 
in the cladding and the coolant. However, in this analysis 
it is assumed that the total energy release is in the fuel 
region. This should not create any significant error. Us- 
ing this assumption, the unsteady heat ection equation 
for the cladding becomes 

oT 


: C 
Ve(k VTE) a Cnc BET 


3. Coolant Region 


Once again, conservation of energy will lead to 


the heat conduction equation plus an additional term which 


a3 





takes into consideration the coolant flow. The governing 
equation is 


( OM ot Creare) as) = 


Ve(keo¥T eg) - Yeo Bz!" co pco co EO CO co ot 


CO co 


where Ve - coolant flow velocity 


Equations (24), (25), and (26) are the governing 

equations used to model the energy transport problem. 
4. Interface Conditions 

The interface between the fuel and the cladding may 
be an actual gap with a finite distance, or the surfaces 
May be in intermittent contact on a microscopic scale. To 
model the heat transfer across this interface, a gap heat 
transfer coefficient 1S introduced. The gap coefficient 
must take into consideration many items (e.g., radiation 
heat transfer across the gap, heat transfer by solid-to-solid 
contact, heat conduction across a gas filled gap). The pre- 
diction of this gap coefficient is extremely complicated and 
beyond the scope of this work. In Ref. [6], a set of values 
oir Hoan is given and the axial variation of Hoan In ies 
analysis is approximately the same. A cosine curve has been 
fitted to the sample data to determine the gap coefficient, 
See Figure 3. The heat flux across the fuel-clad interface 


1s, then, 
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The heat conducted out of the fuel is governed by Fourier's 
equation and is equal to the heat transferred across the gap 
aT, 
q = SS i (28) 


Equating equations (27) and (28) gives the fuel-clad inter- 


face condition 


ke (a,Z) oT 


: : 
Te(a,z,t) + tap) ee = T.(a,z,t) (29) 
apaetrom continuity 
oT. ali 
ke (a,Z) ay—(a,z,t) a kK (a,2) ap (aZ,t) (29a) 


As is common practice in convection heat transfer 
analysis, a coolant surface heat transfer coefficient, hsurf, 
is used to account for the thermal resistance at the clad- 
coolant interface. Similar to the fuel-clad interface analy- 


sis, the clad-coolant interface conditions may be determined 





Kk (b,2Z) ae . | . ) 
T.(b,z,t) + ~—o— ~——(b,z,t) = T Oper ae 30) 
C ane or CO 
and 
aT. alee 
kK .(bsZ) 5p (b.z,t) = Ko (b»Z) ar (baz et) (30a) 


9. Boundary Conditions 


The boundary and initial conditions for the heat 


transport problem are: 
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i) Se (ORs 2) 
yee (r,- Seve T 
Co > (25 PLENUM 
oT 
co H _ 
3) a Cae xt) = 0 
ole 
4) ar Geez cc) = 0 
oT oT 
F H a C H x 
5) Som Chae 5 t) Coa Gatns x» t) =a 
6 ) a es Wels) 
7) ras oe) es ie, 
8) T.9(f.0) = T2 (r) 


Boundary condition 1) results from the assumed azimuthal sym- 
metry. The coolant has been assumed to enter the flow 
channel at a constant temperature, ToL ENUM: This results in 
condition 2). Boundary conditions 3) and 5) result from an 
assumption that no heat is transferred axially from the fuel 
rod. Boundary condition 4) is the result of the use of the 
equivalent cell. Conditions 6), 7), and 8) apa Pac assumed 
Social conditions. 

Solution of the nonlinear coupled neutronic and 
energy transport problem involves the solution of the partial 
Mimnerential equations (11), (12), (21), (24), (25), and 
(26), with the appropriate boundary, interface, and initial 
Sonditions. 

Several works, Refs. [412 311], fice naverdenon— 
Strated the feasibility and success of the finite element 


method in solving nuclear reactor dynamics problems. The FEM 


7 


is used to reduce the partial differential equations developed 
in this analysis to a system of ordinary differential equa- 
tions. Integration of these ordinary differential equations 


MODE) yields the solution. 
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Pie ee CEMENT RORMU AT LON 


In this section, the basic theory underlying the finite 
element method is formulated. Selection of the finite ele- 
ments and the shape functions for the elements are given. 
Some simple transformations which facilitate the integration 


necessary in the FEM are also presented. 


hieeeensolC THEORY 

To obtain a numerical solution, the governing partial dif- 
ferential field equations are transformed into a system of 
ODE in finite dimensional vector space. This may be accom- 
plished in several manners such as the finite-difference 
method, the variational method, or the weighted residual 
method. In this work, the Galerkin method (a weighted re- 
Sidual method) is utilized for the discretization of the 
Spatial domain. The Galerkin procedure will be applied to 
each of the governing field equations. Any of these equa- 
tions may be considered to be in the following form 


se (r,t) - f£y(rst) = f (r,t) (31) 


where w represents the unknown function, e.g., in equation 
(21), w represents ,- , ~ represents the operator for 
Saeh individual equation, and f 1318 a forcing function. In 
the finite element method the solution is approximated as 

N 
tee = eet) = DON (rh) ay 


(t) = <N.>{v43 (32) 
1=] 


1] 
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where N - the number of degrees of freedom 

N. - the element shape function 

Vv. - unknown coordinate function 

< “e-aiatrixapotacivon fora row vec cor 

{ } - matrix notation for a column vector 


Ss - <N, N 


] Dre ey 1] eee N 


Ny 
N 


2 
{N.} = 


de 


The residual, R(r,t), is a measure of the error in this 
finite element approximation. The residual may be considered 


as 
_ ow 
R = AE: a Lv - f (33) 


The best solution for is one which "minimizes" this 
residual. Various "minimums" are obtained by the weighted 


residual method by setting 


fue R dV = 0 Ve big Zack (34) 
V W. (1) - weighting functions 

With the Galerkin method, the weighting functions are the 
Shape functions defining the approximation of equation (32) 
(i.e., We = N.). A noteworthy attribute of the Galerkin 
method is the opportunity of using an integration-by-parts 


40 


of the terms involving the second order spatial derivatives. 
A lower order finite element may be used than would have 
been possible otherwise. Once the weighting functions have 


been chosen, the problem becomes 


f N. (2 - py - f} dV =0  =1,2,...N (35) 
V 


The integration involved in equation (35) is carried out 
on the element level, taking advantage of the use of a "local" 
coordinate system. Once the integration is accomplished, 
the results are merged into a system using "global" coordi- 


nates. On the element level 
ye = <N>8 haa ots) Oe ee me (36) 


where the superscript e indicates the element level. Sub- 
Seating d© into equation (35) and noting {y;}° is not 


a function of the spatial domain yields 


e RR or b e ee é e_ 
f<Ni> (N,}o dv" tw 5} -f [<Nj> £I{N;} -<N > f° ]dV"{p,}° = 0 
V V 

(ee) 
where Tas = 1 oeeeeeans 
N° - number of degrees of freedom for 
an element 
The operator Lo owill vary depending upon which governing 


equation is under consideration. 
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B. SHAPE FUNCTIONS 


The shape functions, N are chosen to satisfy certain 


9 
completeness and convergence criteria [13] and will depend 
upon the finite-element used for the spatial discretization. 
Many previous works, Refs. [4], [11], and [12], utilized 
linear triangular shaped elements to discretize the spatial 
domain. This element was the first element considered. 
However, because the width of the cladding is very thin, 
elements in the cladding region would have extremely large 
aspect ratios (ratio of base to height) unless an extremely 
large number of elements in the axial direction were used. 
A large number of elements becomes numerically untractable. 
Previous experience with triangular elements had shown that 
large aspect ratios yield inaccurate results. Zlamal, Ref. 
[14], showed the error, e, when using triangular elements, 
17S proportional to the square of the longest side, h, and 


inversely proportional to the sine of the smallest angle, y 


e a n°/siny 


A triangular element with a large aspect ratio necessarily 
must have a small related angle which adversely affects the 
error in the FEM. Hopefully to alleviate the problem, an 

1Soparametric, quadratic, rectangular element was selected. 
The aspect ratio would still be large, but experience, Ref. 


[15], with the use of rectangular elements indicated that a 


large aspect ratio is not always a detrimental factor. 
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The shape functions for this element are well documented, 
Ref. [13]. Utilizing a "local" coordinate system (See Figure 


4), the shape functions may be written as 


Corner nodes 1 S="elee oe, 
Ne = (1480) (1400) (Egtno-1) (38) 
Midside nodes N, = 5(1-6°)(14n.) i=2,6 (38a) 
Ne = g(1tE,)(1-n®) i=4,8 (38b) 
where So = aoe 
No a AN 


These normalized shape functions are shown in Figure 5. 
The local and global coordinates are related by the 


following 


5 = <N.>° (ede (39) 


and Z = <Ni>° con (39a) 


C. COORDINATE TRANSFORMATIONS 
When using a local coordinate system, some simple trans- 
formations facilitate the integrations required by equation 


(37). In cylindrical coordinates with azimuthal symmetry 
dV = 2urdrdz (40) 


The derivative terms may be transformed by the following 


ON . ON. 
7 7 
Che ~1) 95 
ne (e CT on, (41) 
Cale on. 
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Figure 4, Element Transformation 
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Figure 5. Normalized Shape Functions 
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where r2ipee is the inverse of the 2x2 Jacobian matrix de- 
fined in Appendix A. As shown in Appendix A, this inverse 


can be easily shown to be (for this problem) 


* * 
“112 
ly}! = (42) 
* * 
Jo, 422 
where 
* * 
Jay = 327 rs Jao = 0 
* = * i 
Joy = 0 ; Joo == 2/72 
Ar - radial length of the element 
QZ - axial length of the element 


Elements of area transform as 
drdz = det[J] d&dn as 


For this particular problem, det[J] may be shown to be 


(Appendix A) 
ne 
det[J] a tg (44) 


where aA© - area of the element 


Elements of axial length become 
an 
dz = om (45) 


where ead a] length of the element 
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Utilization of these transformations makes the integrations 
required in equation (37) amenable to integration by numeri- 


cal Gaussian quadrature. 
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Veen eC AT LONGO R CREite GeeGOn Enanienic 
FIELD EQUATIONS 
In this chapter, equation (37) is applied to the govern- 
ing field equations previously derived. The element matrices 
for each of the operators are developed so that the discreti- 
zation of the spatial domain may be accomplished. The inte- 
grations required by the application of equation (37) are 


performed numerically using Gaussian quadrature. 


A. GAUSSIAN QUADRATURE 
Prior to the application of equation (37), it is appropri- 
ate to discuss briefly the procedure used for the numerical 


integration. The product Gaussian quadrature formula is [16] 


m m 
JW. ie 4 
i (E,n)dndé = x S WW, g(;on,) (46) 
-] -1 =] j=l 
where I, - area integration 


a(&é,n) - any function of & and n 


W. - weight associated with location 
ae ior j 
m - number of Gauss sampling points 


in one-dimension 


The values of the weights associated with each Gauss point 
are given in Ref.. [16]. Equation (42) may be simplified 


somewhat by combining the summations and weights 


I, = We g(E, ony) (47) 
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where k a a 


We - W. X Ws 


For line integrations, the Gaussian quadrature formula in- 


volves only one summation 


m 


| 
1 =f f(n)dn =O Wy F(nj) (48) 
-| 


i=1 


B. NEUTRONIC FIELD EQUATIONS 
The discretization of the spatial domain by the finite 
element method is accomplished by applying equation (37) to 
the governing field equations, using 
= 


d, = <NE> {Hy GTF jis 2 aawemte (49) 


k=F,c,co 


1. Fuel Region 
The governing equation for the fuel region, equation 


(20) after applying equation (37) becomes 


N. 36 a¢ 
1 F 1 @ F 
anf fF yen rardz- 20 f fNi oe spl rOp a) 
i eZ ous 
9 
3 F ° ° 
az (De ae-) 2 poe (1-BI [ko +o-bin(T./Te) | 


t 
—_ = = = y = = 
+ tar. fe AGREE) (Koo)g edt ticee "i rdrdz = 0 
0 


(50) 
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The second order terms in equation (50) may be re- 
duced to a first order term by application of Green's Theorem 
or equivalently integration-by-parts (See Appendix B). 
Dividing through by 2m and reducing the second order terms 


yields 


a 

5 2 7 iE 

J eno, <= az + frnyd— so ar of f A se raras 
Z 


Z if i 


ON. dO¢ ON. do, 
1 z 7 iF ° 
of }PeCar ort oz az TNF ape p (1-8) [kote-bIn(TE/Fe) | 


== SX Ct=t jeans , spo At 
-Wi3Bra, fe (koto )o-dt -N,AC%e rdirdz Os CSae 


From continuity and boundary conditions the line in- 
tegrals are zero. Now using the approximate functions of 
equation (49) and noting that {ye} iS nNOtwa TUNCEITCH Tor 


Space, equation (51) may be written as 
1 ff, ey .@ ef e e e e e 
Vp {N. } alli? rdrdz{v_} ue DELIN; 3 Ne ee aU ot oN a Jrdrdz{p_} 


t 
Mee (1-8) [kotp-bIn(T_/TE) JEN, }°<N.>“rdrdz {pp} -ABE arf ° “Mt-t") (po4p)at! 


@ s,5 -At a 
His, r-rprerdztog - \C°e Af (n,}rare = 0 (52) 
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The integrations may not be easily carried out with 
a shift to local coordinates and with use of the previously 
derived transformations. Rearranging equation (52) and 
assuming the properties are constant for each time step gives 
1 1 


wey f (CIN. 
F By] 1 


* 
N. ~>+{N 


sy Js 


* 
en 


* * 
e 
i n922 }<Jo5 N. pdrdet[d Idédn{p,) 


oe 
Egg (1-8) (Kot0) LF {NG }<N > rdet[J]d&dn{ve}° 


| er 
F e 
ap Bib A In Te th HN pret [ Jdzdn {be} 


ae -, | | 
~\BZ se fac) ss f (N; }<N.>rdet[J]}dedn{p}°-2C%e" ** fel {N. trdet[J Jd&dn 


11 
+ < ff {N.}<N.>rdet[J]dédn {h-}° = 0 (sen 
a J 


Since the element chosen has eight degrees of freedom 
(nodal points), the discretized matrices which result from 
the integration of equation (53) will be 8x8 matrices and the 
forcing function will be an 8xl vector at the element level. 


Defining the matrices as 


1 1 * * * * 
a 7 [{IN; 291, }<Jq, N, >ttNe 1J50 }<Jno N, > Irdet[d] d&dn 
2 

} SAS et «9 
SLHT2; Igyg = a7 De LUNG ed CNS py day 

in 

rn * 
Ne KONG nk Y22 J MK he (oe) 


5a: 








8 
where W. = WW, and ry = )y 
nl ae m2 
e {N. j}<Nyordet[d ]dédn = (H3 «Ig Vg = iE ls (Na) (Ne) rE Wy (55) 
k=] 
Z 
(aaa Os 
! {N, }rdet[J]dédn = {F.}p.. =F 21 (N erly (56) 
ie 
tT p/Tp) tN. }<N. jrrdet lJ ]dédn = (H4 : - Jexg 
2 
ae e M 
=5 AeinTe/Tee) KONG) ONG) PEW (57) 


To carry out the summation of equation (57), the 
temperature Te must be known; however, the temperature 1s 
exactly what 1S being sought. To alleviate this problem, a 
linearization is used. 

In the solution technique, the temperature is pre- 
dicted for the next time step. It is this temperature which 
1s used for the determination of matrix H4. 


Equation (53) simplifies to 


(Dp(H12, «]-[25¢ (1-8) (koto) +482, pF (t) ][H3; ,]+2¢(1-8)b(H4, J} {yp} 


~no At e, ] ae 
=e 1C {Fs} 2 ylh3 5 {be} 0 C539) 
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The function f(t) is evaluated by summing the 
Values at each time step using the trapezoid rule for numeri- 


cal integration. 


ree | Lene 
f(t) = e7** 5 a7 *® Fxeto]de'! (59) 
0 
Defining 
intl 
I.{g(t) ] i 2 h.{g(t.)+ Gita ae (60) 
7Ct) e (koto) 
h. - time step taken 
din ae 


The function may be expressed as 
¢ | 
_ _-At ; 
f(t) = e 2X Le igem (61) 
T= 


S 


number of time steps 


2. Clad Region 


Following the same procedure as with the fuel equa- 


mon. the discretized form of equation (11) becomes 
e, | * ,e _ 
{0 fH12; 5] i Za LH3; ,J}{,} i Weave +0 (62) 


3. Coolant Region 


The governing equation for the coolant region, equa- 


tion (12), may be discretized into the form 


e | : ez 
{D. fH12; 5] - Bacon ty li cg! + pence vee = 0 (63) 
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Once the governing equations have been discretized 
at the element level, they are combined into a system of 
equations at the global level. On the global level the 
governing equation for neutron transport takes the form of 


ie lng] tobe Ube * lite On ea eoee 


nx 
where [H], [P], and [F] represent the system matrices and 
n 1S the number of nodal points used in the discretization. 
There are, then, n simultaneous ordinary differential equa- 


tions used to describe the neutron transport problem. 


Co. HEAT TRANSPORT FIELD EQUATIONS 

The spatial domain for the heat transport problem is 
discretized in the same manner as the domain for the neutronics 
problem. The same element matrices previously defined are 


Valid. Let 


. e e aes 
ie . se tae J ligZaeeen (65) 


K=P 565 CO 


me Fuel Region 


The governing field equation for the fuel, equation 
(24), is discretized by applying equation (37). Using an 
integration by part to lower the order of the second order 
terms allows equation (24) to be written as 


oT 


eae a oreo Se ae 
E [N= k. Sala ote [N.k- oF ae a {kel or 1 acre. 
oT. 
“eZ ecN. b- “- Ni0pC oe ae rdrdz = 0 (66) 
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From continuity considerations the line integrals 
are zero except along the boundaries of a region. It is 
assumed that no heat is transferred from the cell in the 
axial direction (boundary condition 5); therefore, the 
first line integral of equation (66) is zero. In the neu- 
tron diffusion problem there was continuity of flux at the 
interfaces so that the line integrals were zero; however, 
the heat transfer at the interfaces is affected by the gap 
and film conductances. The fuel-clad interface condition, 
equation (29), may be rewritten as 
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Substituting equation (67) into (66), dividing by minus one 
and utilizing equation (65) yields 
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Transforming to local coordinates and integrating 
by Gaussian quadrature yields the same element. matrices as 
given for the neutron flux, equations (54) and (55), except 
for the line integrals between regions. It should be noted 
that the line integrals exist only on the interfaces, along 
which there is a discontinuity of temperature. For the 
fuel equation the interface corresponds to the local co- 
ordinate €=1]1. Define 


2 


1 e 
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where the N's are evaluated at —=1. Many of the terms of 
Kl will be zero since only the nodes on the &=1 boundary 
will have shape functions which are non-zero. It is through 


Kl that the temperatures for each region are coupled to- 


gether. Equation (68) may now be written as 
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~eDee [H3]{yp}° + opC - [HB] {t,} (70) 
In obtaining this equation, it was assumed that material pro- 
perties for each nodal point were constant at each time step. 
Perhaps a better assumption would have been to assume an 
average value for the properties of each element. The dif- 
ference should not be significant, and the assumed constant 


nodal properties were numerically more tractable. 
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ee Glad Region 
Applying equation (37) to the governing field equa- 
tion for the clad, equation (25) gives the discretized form 
of the equation. Assuming no heat transfer in the axial 
direction on the boundaries (Boundary condition 5), equation 


(25) becomes 


3T ANE oe Ale aE 
Gib 1 Cc 1 €! 
~ FUNK. Bplay dz + SS klar a + ae az] vardz 
Z az 
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ee. Pe aye N. st rdrdz = ) (71) 


In the cladding, there are two interfaces along which 
the line integral of equation (71) is not zero, along the 
fuel-clad interface and along the clad-coolant interface. 

For the fuel-clad interface, equation (61) applies. For the 
clad-coolant interface, the interface condition, equation 
(30), may be rewritten as 
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Along the fuel-clad interface, the local coordinate 
corresponds to €=-1. Define the new element matrix 
i | : 2 fem | 
A f {N,}°<N,>"dn = [kK2; slog = 3 de (Nady (Na), Wy 
7 k=] 
where the N's are evaluated at §=-1. As with Kl, Ke will 
have many zero values because the shape functions are 


evaluated at &=-1. 
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Along the clad-coolant interface, the local coordi- 
Meee corresponds to €=1 and the K1 matrix is appropriate. 


Substituting equations (67) and (72), equation (71) 


becomes 
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The governing equation for the clad region may now be written 


as 
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The line integrals of equation (73) affect only nodes which 
are on one of the boundaries; therefore, the nodal inputs 
into Kl] and K2 are zero unless the node is on one of the 
Doundaries. 
3. Coolant Region 
The field equation governing the coolant may be 
discretized in the same manner as above. After applying 


the Galerkin method and performing an integration Dy parts 
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on the second order terms, equation (25) 


becomes 
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The line integral, when evaluated at c, is zero (boundary 
condition 4). When evaluated at b, or correspondingly at 
meal, equation (72) is valid and K2 matrix is appropriate. 
All the terms of equation (75) have been defined except the 
flow term. Define 

in 
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Transforming to local coordinates and integrating reduces 


Baiation (75) to 
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Now that the governing equations have been defined 
for each region on the element level, equations (70), (74), 
and (77), they may be assembled into a system equation on 
the global level. The equation will be in the general form 


of 
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Wee DISCRETIZATION OF THE SPATIAL DOMAIN 

Prior to the numerical solution of the governing equa- 
tions, equations (64) and (78), the spatial domain must be 
divided into a number of elements. For this work the domain 
was subdivided as shown in figure 5. 

Since there is a discontinuity of temperatures at the 
interfaces, as described by equations (29) and (30), a 
noval application of the FEM method was necessary. The com- 
mon practice for handling these "flux" type boundary condi- 
moms 1S to define a constant reference temperature, 7, ; 
as when working with a convection heat transfer problem [17], 
or to define a Known function, as when working with a 
fracture mechanics problem [18]. In either case the refer- 
ence condition was known. The novel application here lies 
in the use of a different field equation to describe the 
reference temperature, e.g., the clad equation (74) is the 
reference condition for heat transfer from the fuel across 


the gap interface. 
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The discontinuity of temperature at the interface neces- 
Sitated another novel application of the FEM. Since there 
1S a temperature drop along each interface, a single node 
there 1s not adequate. In the discretization of the domain, 
two nodes were used for each interface point (for example, 
Points 62 and 63 in figure 6). This allows the temperature 
drop due to the gap and film conductances to be taken into 
consideration. Since two nodes are used, the governing equa- 
tions for each region are not directly coupled together. 

The coupling of the regions is accomplished by the "flux" 
boundary or interface conditions since it 1s assumed that 
any heat flux leaving a region enters the adjacent region. 


Consider a typical set of elements on an interface 
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Figure 6. Finite Element Discretization 
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The coupling terms Kl and K2 may be combined into a system 
K matrix which shows the coupling. The K matrix for the 


Simple set shown is 
pu 4a 27627) 8 oe 2913 14-15. 16 
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As can be seen, the nodes on the interface are coupled to 
the adjacent element interface nodes. For example, node 2 


in element I is coupled to nodes 3, 8, and 14 in element II. 


meer riMUM COMPACTING SCHEME 

The system matrices (H, P , etc.) are nxn matrices, 
where n is the number of nodal points used in the discre- 
tization of the domain. In terms of computer storage, these 
matrices may become excessively large if they are stored as 
nxn. There are several techniques available to reduce this 
Storage requirement. The most common method is the banded 
Storage scheme, whereby only the banded portion of the 


matrices are stored. With judicious numbering of the nodes, 
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considerable savings may be realized. However, it is not 
the optimum storage scheme [8]. 


th 


Since the shape functions, N for the k—— nodal equa- 


ce 
tion are nonzero over only the element containing k, the 
System matrices are not only banded but sparse as well. 

The sparseness is due to the non-consecutive numbering of 
the nodes surrounding the yall node. The optimum compact 
storage (0CS) scheme compacts the matrices by storing only 
the non-zero elements of the matrices. The implementation 
of the OCS scheme requires two additional integer arrays, 
say JA and NAME. The NAME array identifies the nodal points 
which contribute to each nodal equation. The JA array acts 
as a pointer to indicate where the nodal equation starts in 


NAME. Consider the following simple 2x2 system with nodes 


as indicated 





The NAME array starts with node |] and identifies the nodes 
which contribute to node 1] (i.e., 5, 4, and 2). The NAME 
array would, then, give the nodes contributing to node 2 
ema so forth, so that 
Woo e656 789010 + 1i 1219 14+ «.. % JC 
auetieemeec5a2 > 254163 ., 365 2 2 ss. « 9656> 


and ‘meee 
ape ae =<] 5 11>... > 
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The algorithm to assemble the element matrices into a compact 
storage vector is straightforward and represents a significant 
Savings in computer storage [8]. The system matrices are 
Stored as a vector rather than a two-dimensional array. For 
example, the value which would be stored in position (1,5) 


of the nxn array is stored in position 2 of the system vector. 
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VI. NUMERICAL SOLUTION 


This section contains a brief description of possible 
solution techniques in addition to the solution technique 
chosen. Computer subroutines necessary to implement the 


technique are also described. 


fee ELECTION OF METHOD 

The numerical solution of the system of implicit ordi- 
nary differential equations, equations (64) and (78), may 
be accomplished by any of a number of different techniques 
such as Houbolt's method, Crank-Nicolson's method, Gear's 
method, or implicit Gear's method. It was not the objective 
of this analysis to determine which of the numerical solu- 
tion schemes is the most efficient. Each method has its 
advantages and disadvantages. The Crank-Nicolson method 
is a single-step, implicit equation solver and, therefore, 
does not require storage of previous time solutions. When 
analyzing neutronic problems, the system of equations which 
arises is commonly very stiff (i.e., a rapid change in flux 
Over a short period of time). The Crank-Nicolson method 
has, in a past work [8], demonstrated difficulty in tracking 
these stiff systems. Gear's method was specifically 
developed for stiff systems and can handle the problem very 
well. However, Gear's method is a multi-step, predictor- 
corrector method requiring storage of previous time solutions. 


In addition to this disadvantage, Gear's method requires the 


66 





transformation of the developed implicit 0.D.E.'s into an 
explicit system of 0.D.E.'s. After this transformation is 
done, the system matrices are no longer sparse or banded, 
thus eliminating the use of the optimum compacting scheme. 
In an effort to overcome these difficulties, Gear's method 
was modified, Ref. [9], to treat the implicit system of 
equations as well as to allow use of the optimum compacting 
Scheme. A previous work, Ref. [8], has Shown that the 
implicit Gear's method is particularly attractive in solv- 
ing the type problem developed in this analysis. Therefore, 
the implicit Gear's method is used for the solution of the 
System of 0.D.E.'s arising inthis analysis. 

No attempt will be made here to give the mathematics in- 
volved in developing the implicit Gear's method. Reference 
[9] may be consulted if details are desired. A listing of 
the computer program developed will be given in the Computer 
Program section. In order to utilize the implicit Gear's 
method, several user supplied subroutines must be developed: 
1) DIFFUN, 2) JACMAT, and 3) NUITSL. 

B. USER SUPPLIED SUBROUTINES TO IMPLEMENT 
fae IMPLICIT GEAR*S METHOD 
om, OILFEFUN 
Subroutine DIFFUN evaluates equations (64) and (78) 
for a given time and for given values of wv, w, t and Tt 
Bimece at each nodal point, i, there is a solution for the 
flux and for the temperature, the solution was set equal to 


DYI and DYII, respectively. In addition to having flux and 
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temperature at each nodal point, there are also three differ- 
ent regions in the domain which have different governing 
equations. An-integer array, ITYPE, was developed to indi- 
cate for each nodal point whether it was: 0) a fuel node 

not in an interface element, 1) a fuel node in an interface 
element, 2) a cladding node or, 3) a coolant node. Using 
ITYPE, the computer program is directed to a different sec- 
tion depending upon the type of node being considered. After 
all the nodes have been considered, boundary conditions are 
established by changing DYI and DYII for the appropriate 
boundary nodes. Since in this analysis there is continuity 
of flux at the interfaces, special considerations must be 
Given to these nodes. At the fuel-clad interface, the value 
of DYI for the clad node was set to the value of the flux 

at that node minus the value of the flux at the adjacent 


- yp. 


node (i.e., DYI, = a2 


1): Similarly, at the clad-cool- 
ant. interface, the value of DYI for the coolant node was 

set to the value of the flux at that node minus the value 

of the flux at the adjacent node. During the solution of 
the problem, DYI is driven toward zero, which in the limit 


forces v; to equal wy This is the desired continuity 


i-] 
Result. 
2. JACMAT 
Subroutine JACMAT, evaluates the Jacobian matrix 
(for Gear's method) at the given time and for the current 
values of the dependent variables. The Jacobian for an equa- 


tion of the type, 
vey ste ae (79) 
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may be represented as [19], 
(80) 


where Ot 6 and By are coefficients from Gear's method and 


h is the time step. Using the notation of DIFFUN, let DYI 
and DYII represent equations (64) and (78), respectively. 
The Jacobian matrix may, then, be written as (J is called 


PW in JACMAT.) 


elt | so sD envi | se 9DYIT, 
OW Bh eo Boh 


—= (81) 
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It is the form of equation (81) which is programmed in 
JACMAT. As in DIFFUN, the integer array ITYPE is used to 
indicate the appropriate section of the program to be 
utilized. The problem boundary conditions must also be 
accounted for in JACMAT. In DIFFUN, the value of DYI or 
DYII was set to zero for constant boundary conditions (i.e., 
zero). This cannot be done in JACMAT since a division by 
zero would occur. For a constant boundary conditton at the 
ith node, the value of PW is set to one for the diagonal 
term and zero for all other terms of the ith equation. 
See NUL TSL 

Subroutine NUITSL solves the system of equations 
for the quasi-Newton iterates. In this analysis the system 
is solved using a successive over-relaxation (SOR) method. 
In this work, the optimum amount of over-relaxation was not 


determined. Since no effort was made to find the optimum, 
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it was felt a small over-relaxation would be best. The over- 
relaxation factor of 0.02 was used. For small values of 
this factor, the SOR method approaches the Gauss-Siedel 


iteration technique. 
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VIT. PROCEDURE 


In this section, the method utilized to obtain a solu- 
tion is described. The input data necessary to run the 
developed computer program will be documented. 

Prior to initiating a transient overpower excursion, the 
Steady-state conditions for the fuel cell must be known. 
Since the systemof equations which were developed are not 
Specifically designed to obtain a steady-state solution, the 
initial steady-state conditions must be part of the input 
data. The initial temperature distribution was obtained 
from the steady-state conditions given in Ref. [7]. 

The axial temperature distribution for the fuel center- 
line, fuel surface, clad, and coolant have been determined 
for several different fuel life cycles [7]. For this analy- 
Sis, the beginning of life cycle for channel 10 was used. 
eeenough this distribution is somewhat artificial, it should 
be adequate for this analysis. It is the trends of the 
results which are considered important. The distribution 
within the fuel radially is taken to vary as the square of 


the radial distance; then 


re re 
Te(r,z,0) T.(0,2,0)(1- wo) - Tpla,z,0) (3) 


Within the cladding and the coolant, the initial radial tem- 


Perature distribution is assumed to be constant. 
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The initial flux distribution is assumed to be radially 
constant, a flat flux assumption. In the axial direction 
the flux is assumed to vary as the shape of the sine func- 
aon. the maximum flux, the flux at the axial center, is 
an input parameter. For this analysis, the maximum initial 
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flux was taken to be 10 
To obtain a steady state flux distribution, the value 
of fission cross section for the fuel is varied. A trial- 
and-error method is used until a critical fission cross sec- 
tion, ie » which gives a steady flux is obtained. 
Once the steady-state conditions have been determined, 


the excess reactivity may be inserted. This starts the 


transient overpower excursion. 


A. INPUT DATA 

The first data card contains: the order of Gauss quadra- 
ture, the number of radial elements in the fuel, the number 
of axial elements, number of nodal points in the radial 
direction, and the height of the fuel rod. The next cards, 
One for each radial nodal point, contain the nodal radial 
distances. The next cards contain the fuel centerline, fuel 
Surface, clad and coolant temperatures. There is one card 
For each axial node. The next card contains the maximum 
Flux. The next four cards contain the physical parameters 
listed in Table I. The last input cards contain the time, 
end time, estimated initial time step, minimum time step, 


and maximum time step. A sample data deck is shown in figure 
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TABLE I. 


Fuel Diffusion Coefficient (DCF) 
Doppler constant (B) 

Energy release per fission (E) 
Fraction of delayed neutrons (BETA) 


Fraction of delayed neutrons for 
the ith group (BETAI) 


Decay constant for ith delayed neutron 
group (DCLAMI) 


Initial flux for delayed neutrons 


Average number of neutrons released 
per fission (ANU) 


Neutron velocity (VEL) 
Fuel absorption cross section (SIGAF) 


Critical fuel fission cross section 
(SIGFF) 


Blanket absorption cross section (SIGAB) 
Blanket fission cross section (SIGFB) 
Step reactivity input (RHOA) 

Ramp reactivity input (RHOB) 

Fuel density (DENF) 

Clad diffusion coefficient (DCC) 

Clad specific heat (CPC) 

Clad density (DENC) 

Clad thermal conductivity (TKC) 

Clad absorption cross section (SIGAC) 
Coolant diffusion coefficient (DCCO) 
Coolant absorption cross section (SIGACO) 
Coolant flow velocity (VCO) 

surface heat transfer coefficient (HSURF) 


ae 


Physical Parameters 


0.93 [cm) 
0.006 
Werle (th 
0.0064 


~l2 [cal/fissions ] 


0.0064 


0.0/84 
1x10 9 [neutron/cm¢sec] 


2.44 

4.8 x10° [cm/sec ] 
0.088 [em™!] 
0.0586875 [cm™ |] 


0.0 rom™ |} 

0.0 fom |] 
Variable 

Variable 

10.9 [gm/cm? ] 
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0 [gm/cm?] 

.0526 [cal/cm sec °C] 
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2S [eh 

.00004 Com |] 
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0.7 [cal/em¢sec °C] 
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Vee cles 


When using the finite element method, one of the first 
considerations must be given to the convergence of the 
method. To determine convergence, the results for a given 
point are compared for different finite element discretiza- 
tions. As shown in figure 8, the results are comparable 
but no definite claim of convergence can be made. However, 
for this work, it was felt that these results were adequate. 
It was not the object of this analysis to arrive at the 
"final" result; it was the trends and methods that were of 
interest. Since the 66-element mesh appears to give a fair 
approximation of the results, the 66-element mesh was used 
as the discretized domain. 

The next item of consideration was the determination of 
a neutronic steady-state condition. This proved to be a 
very time consuming task. The fission cross section for the 
fuel was varied by a trial-and-error method in an attempt 
to find the critical cross section which would give a steady 
State. As may be seen in figure 9, a change in cross sec- 
tion of less than one percent significantly affected the 
State of the problem. It was felt that the critical value 
was between the values of 0.05875 and 0.058625. Time did 
not permit investigation for critical value; therefore, 
it was assumed that the value for the critical fission cross 

er 


section was half way between the values (i.e., Le = GeUsee e730). 
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Even if this value is in error, which it most probably is, 
the net effect would only be a small decrease or increase 
in the proposed reactivity insertion. The reactivity inser- 
tion would, then, be an approximation of the actual reactiv- 
ity of the problem. 

The first test problem considered was a step increase 
in reactivity of approximately ten dollars. For the uranium 
dioxide fuel, one dollar of reactivity was taken to be 
0.0064. Figure 10 shows the time history of the flux at 
the center of the fuel rod. Figure 11 gives the correspond- 
ing temperature profile. The temperatures were taken at the 
hottest point of each region at the axial center (i.e., the 
fuel centerline, clad inside surface, and coolant inside 
surface). As seen on the fuel temperature time history, the 
fuel rapidly reaches the fuel melting point. The model 
developed does not take into consideration melting of the 
fuel. This melting would tend to decrease the effect of the 
transient. The problem was allowed to continue despite this 
inconsistency in the mathematical model. A short time after 
fuel melting, the inside surface of the clad reaches its 
melting point. The temperature in the coolant experienced 
What is felt to be a numerical phenomenon. The coolant tem- 
Perature decreased prior to the small rise at the end of the 
transient. Intuitively, this decrease does not seem to be 
realistic. A similar occurrence was observed while conduct- 
ing sample tests on the developed computer program. Refer- 


ence[20] reported the same phenomenon. It is felt that this 
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phenomenon is a quirk of the finite element method. The 
radial and axial distributions of the neutron flux are pre- 
sented in figures 12 and 13 for time equal to /.39 x nay 
Seconds. The distributions are, basically, as anticipated. 
The neutron flux peaks slightly before the axial center. 

It was expected to peak at the axial center. The radial 

and axial temperature distributions for the same time are 
presented in figures 14 and 15. As with the flux, the tem- 
perature profiles were, basically, as expected. The fuel 
temperature peaks slightly below the expected location, 

most likely in response to the peak in axial fluxes. The 
coolant unexpectedly drops near the outlet of the fuel rod. 
The finite element method characteristically has some prob- 
lems on the boundaries of the domain; this may account for 
the drop in coolant temperature. 

The temperature of the fuel and cladding do not appear 
to be as closely coupled as anticipated. As seen in figure 
14, a significant increase in fuel temperature has resulted 
in a relatively small clad temperature increase. As noted 
in figure 11, there appears to be a time lag in temperature 
response for each region which may account for part of the 
apparent temperature disagreement. It is felt that the 
temperatures should be more closely related, which indicates 
a higher gap heat transfer coefficient should be utilized. 
Since values of the gap heat transfer coefficient were 
assumed, it is not unreasonable to believe the values used 


are too low. 
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Time did not permit investigation of other reactivity 
insertions. Other reactivity inputs may be investigated by 
Students in the future. 

This work does not represent a solution to the very com- 
plicated nuclear reactor problem. It does represent an 
application of a numerical technique which is relatively new 
to nuclear applications. Methods for implementing the finite 
element method have been discussed, and a computer code has 


been developed for the simplistic model considered. 
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IX. RECOMMENDATIONS 


For the model developed, perhaps the most important item 
tO pursue iS the critical fission cross section. A better 
determination of this value is necessary so that the reactiv- 
ity insertion is more accurately known. Different test cases 
Mometne prompt critical and prompt subcritical reactor could 
then be conducted. 

In further developing the model, more consideration 
Should be given the gap heat transfer coefficient. As noted 
in the results, the value used appears to be too small. 
Sample problems for different gap heat transfer coefficients 
would give a better indication of the values to uSe. 

Melting of the fuel during the transient would probably 
be the next major improvement on the model. With relatively 
few changes, the model could be adapted to allow melting 
element by element. This, too, would be an approximation 
but, still, an improvement to the model. Perhaps at the 
Same time, a simplified model to take into consideration the 
fuel restructuring could be implemented. 

Another improvement would be to consider reactivity 
feedbacks in addition to the Doppler feedback. Sodium void- 
ing and fuel rod expansion are two of the more important 
feedback effects to consider. 

On the numerical side, probably the most important thing 


to do would be to run the computer program on the 


Sy 





"H-compiler", which optimizes the program. However, on 
Several runs using the H-compiler, erroneous results ven 
Obtained. With sufficient time, this could be corrected 

to allow use of the H-compiler. The use of the H-compiler 
results in a savings in computer time. The present program 
runs on a "G-compiler" and takes excessive amounts of com- 
puter time (two to four hours per run). 

In addition to this, the optimum over-relaxation factor 
in the implicit Gear's method could be determined by trial- 
and-error. 

Implementation of these recommendations should enhance 


the analysis and lead to a more efficient computer code. 
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APPENDIX A 
DEVELOPMENT OF TRANSFORMATIONS 


The Jacobian matrix [J] may be written (8) for two 


dimensions as 
N 


i,cct z Ni 2; 


T=] 


[J] = (Al) 
N 


N 


M= 


— 
WN 
<i 


> NG hy > Ne onal 


i= i=] 


For a simple 2x2 matrix [A] the inverse is 


a b 
mf 
C d 


14 d  =b 
LAL © @ettal | _. 


Applying this fact to equation (Al) gives 


N N 
2,8 - 2 Piece 


i=] 1 
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From matrix algebra 
N 
det[J] = DUNG er 
i=] 


The derivatives of the shape functions may be found from 


equations (38). 


] 


Ny og = gl l#n) (2é4n) Nyon 7 gilt) (2nt6) 

Nog = 78(14n) No = p1-8°) 

Ng ¢ 2 g(l4n) (2E-n) Nan = g(1-8) (2n-8) 

Nae 27 p(n?) Nay 2-1-8) (A4) 
Neg = g(T-n) (2E+n) Ne oy = gle) (2nté) 

Nee = 76 (1-n) Neon 7 gI8*) 

Nag = g(I-n) (2E-n) No, = q(18) (2n-6) 

Nee. 5(1-n°) pen a eulice 


If one now considers an arbitrary element 


Zz 









6 
j i | 
| 
| 
| 
{ | | 
r r r 


2 3 r 
with the midside nodes exactly at the midpoint (not a neces- 


Sary criteria for the FEM) and substitutes into equation (A3), 
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mae result 1S 
( A 
detiJ] = —-———"———"—_ = 1. (A5) 


Substituting into equation (A2) and using (A5) will yield 


J ge H(z 26 | Swap le (A6) 
1] ae Deg | So lee 
* 
Jio = QO, (A7) 
* 
Joy ae 5) (A8) 
and is = “2 [s(ra-r,)] = €/Z4-Z, ; (A9) 


The inverse of the Jacobian matrix now becomes 


2/ra-¥r, 0 
tayo! = (A10) 


0 2/Z2-Z, 


For integration along a line, the transformation used 


1s 
dz = det[J' ]dn (All) 
miecnis case 
N 
! — 
ae al? | > NG nti (A12) 
T=] 


Again considering the arbitrary element and substituting 


(A4) into (A12) will result in 


det[J'] = = > (Ads) 














APPENDIX B 
REDUCTION OF SECOND ORDER TERM 


The second order term of the governing field equations 
may be reduced to first order by integration by parts. 


Consider, for example, 


Sf Nex 4 or (rd 2 <b) + 2 (0 oe) rdrdz (B1) 


which may be expanded as 


2 
Sf tw,r0 2s Va 2 ono Me nro AF + nr a drdz 
1 ype 7 dZ 02 


(B2) 


Integrating just the second order terms by parts will yield 
if N.D ary rdrdz =f [N rp oY dz 
17 Re es laws 
or 
Z 
3D ah 
Sf 3 SECON. + rN, = + rd a] drdz  (B3) 


and 


SS D zy rdrdz J p oY) rdr Sf ea N, 22) rardz  (B4) 


mmipStituting the results of (B83) and (B84) into equation (B82) 


Will give 


oN. 
fo rD aries * +P IND oe) dr Sf Diet = + —i3 2] rdrdz_ (B5) 
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APPENDIX C 
LIST OF RELATIONS FOR MATERIAL THERMAL PROPERTIES 


BUEL ( U0.) 


1. Specific Heat, Ref. [19] 
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[cal/gm °C] 


m= 26 
2. Thermal Conductivity, Ref. [5] 
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Properties are assumed to be temperature independent, 
and average values from Ref. [5] were used for the 
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